#relevant packages

rm(list=ls())
require("tidyverse")
require("ggplot2")
require("stargazer")
require("xtable")
require("texreg")
require("reshape2")
require("RItools")
require("estimatr")
require("fastDummies")
require("ggpubr")
require("estimatr")
require("mediation")
require("car")
require("ltm")
require("table1")
require("psy")
require("ggstatsplot")
require("ggplot2")
require("quanteda")
require("readtext")
require("tableone")
require("quanteda.textplots")
require("htmltools")
require("devtools")
install_github("naoki-egami/exr", dependencies = TRUE)
require("exr")
require("kableExtra")
require("patchwork")


### read data
data <- read_csv("data.csv")


## ----------------------------------------------------------------------
# Run main pre-registered models and generate main figures and tables
# ----------------------------------------------------------------------

#Test H1a-H1b (Does shaming decrease support for target?)#

##H1a: favorability
H1a<- lm_robust(target_favorability~shaming, data = data)%>% tidy %>% 
  mutate(.,
         outcome = "Favorability \n (Target)",
         Towards = "target") 

##H1bA: military cooperation
H1bA<- lm_robust(target_military~shaming, data = data)%>% tidy %>% 
  mutate(.,
         outcome = "Military Coop. \n (Target)",
         Towards = "target") 

##H1bB: economic cooperation
H1bB<- lm_robust(target_economic~shaming, data = data)%>% tidy %>% 
  mutate(.,
         outcome = "Economic Coop. \n (Target)",
         Towards = "target") 

#add controls
H1a_controls<- lm(target_favorability~shaming+Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state, data)
summary(H1a_controls)

H1bA_controls<- lm(target_military~shaming+Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state, data)
summary(H1bA_controls)

H1bB_controls<- lm(target_economic~shaming+Democrat+Republican+Independent+
                     Educated+Black+White+Asian+Native+Other+Religious+
                     Female+Age+Political_interest+Conservative+hawk1+hawk2+
                     hawk3+hawk4+state, data)
summary(H1bB_controls)


#Table
#save table of interaction terms
#Create table
tab1<-stargazer(H1a_controls, H1bA_controls, H1bB_controls,
                out= "Tables/H1control.tex",  style="qje",
                title="Testing H1a-H1b with Control variables",
                keep = "shaming",
                dep.var.labels = c("Favorability (Target)", "Military Coop. (Target)", 
                                   "Economic Coop. (Target)"),
                covariate.labels ="Shaming",
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:H1controls",
                # star.char = c("", "", ""),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{4}{c} {\\parbox[t]{14cm}{ \\textit{Notes:}
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex

tab1 <- gsub("\\begin{tabular}"
             ,"\\resizebox{1.1\\textwidth}{!}{\n\\begin{tabular}"
             ,tab1
             ,fixed=T)

tab1 <- gsub("\\end{tabular}"
             ,"\\end{tabular}}"
             ,tab1
             ,fixed=T)

cat(tab1, sep = "\n")

write(tab1, 'Tables/H1control.tex')



#Test H1c (Does the relationship with the shamer moderate H1a, H1bA, or H1bB?)
#create sub datas for Shamer ally and enemy conditions
sham_ally <- data[ which(data$shamer_ally==1), ]
sham_enemy <- data[ which(data$shamer_ally==0), ]

##favorability
H1c_Aally<-lm_robust(target_favorability~shaming,sham_ally)%>% tidy %>% 
  mutate(.,
         outcome = "Favorability \n (Target)",
         Towards = "target",
         Shamer = "U.S. Ally") 

H1c_Aenemy<-lm_robust(target_favorability~shaming,sham_enemy)%>% tidy %>% 
  mutate(.,
         outcome = "Favorability \n (Target)",
         Towards = "target",
         Shamer = "U.S. Enemy") 


##military
H1c_Bally<-lm_robust(target_military~shaming,sham_ally)%>% tidy %>% 
  mutate(.,
         outcome = "Military Coop. \n (Target)",
         Towards = "target",
         Shamer = "U.S. Ally") 

H1c_Benemy<-lm_robust(target_military~shaming,sham_enemy)%>% tidy %>% 
  mutate(.,
         outcome = "Military Coop. \n (Target)",
         Towards = "target",
         Shamer = "U.S. Enemy") 


##economic
H1c_Cally<-lm_robust(target_economic~shaming,sham_ally)%>% tidy %>% 
  mutate(.,
         outcome = "Economic Coop. \n (Target)",
         Towards = "target",
         Shamer = "U.S. Ally") 

H1c_Cenemy<-lm_robust(target_economic~shaming,sham_enemy)%>% tidy %>% 
  mutate(.,
         outcome = "Economic Coop. \n (Target)",
         Towards = "target",
         Shamer = "U.S. Enemy") 

#save interaction models, save in table form later
##H1c_A:favorability
H1c_A<-lm(target_favorability~shaming*shamer_ally,data)

##H1c_B:military cooperation
H1c_B<-lm(target_military~shaming*shamer_ally,data)

##H1c_C:economic cooperation
H1c_C<-lm(target_economic~shaming*shamer_ally,data)


#add controls
H1c_Acontrols<-lm(target_favorability~shaming*shamer_ally+
                    Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state,data)

H1c_Bcontrols<-lm(target_military~shaming*shamer_ally+
                    Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state,data)

H1c_Ccontrols<-lm(target_economic~shaming*shamer_ally+
                    Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state,data)


#Test H1d (Does the relationship with the target moderate H1a, H1bA, or H1bB?)
#create sub datas for ally and enemy conditions
tar_ally <- data[ which(data$target_ally==1), ]
tar_enemy <- data[ which(data$target_ally==0), ]

##favorability
H1d_Aally<-lm_robust(target_favorability~shaming,tar_ally)%>% tidy %>% 
  mutate(.,
         outcome = "Favorability \n (Target)",
         Towards = "target",
         Target = "U.S. Ally") 

H1d_Aenemy<-lm_robust(target_favorability~shaming,tar_enemy)%>% tidy %>% 
  mutate(.,
         outcome = "Favorability \n (Target)",
         Towards = "target",
         Target = "U.S. Enemy") 


##military
H1d_Bally<-lm_robust(target_military~shaming,tar_ally)%>% tidy %>% 
  mutate(.,
         outcome = "Military Coop. \n (Target)",
         Towards = "target",
         Target = "U.S. Ally") 

H1d_Benemy<-lm_robust(target_military~shaming,tar_enemy)%>% tidy %>% 
  mutate(.,
         outcome = "Military Coop. \n (Target)",
         Towards = "target",
         Target = "U.S. Enemy") 


##economic
H1d_Cally<-lm_robust(target_economic~shaming,tar_ally)%>% tidy %>% 
  mutate(.,
         outcome = "Economic Coop. \n (Target)",
         Towards = "target",
         Target = "U.S. Ally")  

H1d_Cenemy<-lm_robust(target_economic~shaming,tar_enemy)%>% tidy %>% 
  mutate(.,
         outcome = "Economic Coop. \n (Target)",
         Towards = "target",
         Target = "U.S. Enemy")  


#report interaction terms of H1d in table
##H1d_A:favorability
H1d_A<-lm(target_favorability~shaming*target_ally,data)

##H1d_B:military cooperation
H1d_B<-lm(target_military~shaming*target_ally,data)

##H1d_C:economic cooperation
H1d_C<-lm(target_economic~shaming*target_ally,data)

#with controls
H1d_Acontrols<-lm(target_favorability~shaming*target_ally+
            Democrat+Republican+Independent+
            Educated+Black+White+Asian+Native+Other+Religious+
            Female+Age+Political_interest+Conservative+hawk1+hawk2+
            hawk3+hawk4+state,data)

H1d_Bcontrols<-lm(target_military~shaming*target_ally+
                    Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state,data)

H1d_Ccontrols<-lm(target_economic~shaming*target_ally+
                    Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state,data)


#Test H2a-H2b (Does shaming increase support for *shamer*?)#

##H2a: favorability
H2a<- lm_robust(shamer_favorability~shaming, data = data)%>% tidy %>% 
  mutate(.,
         outcome = "Favorability \n (Shamer)",
         Towards = "shamer") 

##H2bA: military cooperation
H2bA<- lm_robust(shamer_military~shaming, data = data)%>% tidy %>% 
  mutate(.,
         outcome = "Military Coop. \n (Shamer)",
         Towards = "shamer") 

##H2bB: economic cooperation
H2bB<- lm_robust(shamer_economic~shaming, data = data)%>% tidy %>% 
  mutate(.,
         outcome = "Economic Coop. \n (Shamer)",
         Towards = "shamer") 

#add controls
H2a_controls<- lm(shamer_favorability~shaming+Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state, data)
summary(H2a_controls)

H2bA_controls<- lm(shamer_military~shaming+Democrat+Republican+Independent+
                     Educated+Black+White+Asian+Native+Other+Religious+
                     Female+Age+Political_interest+Conservative+hawk1+hawk2+
                     hawk3+hawk4+state, data)
summary(H2bA_controls)

H2bB_controls<- lm(shamer_economic~shaming+Democrat+Republican+Independent+
                     Educated+Black+White+Asian+Native+Other+Religious+
                     Female+Age+Political_interest+Conservative+hawk1+hawk2+
                     hawk3+hawk4+state, data)
summary(H2bB_controls)


#Table
#save table of interaction terms
#Create table
tab1<-stargazer(H2a_controls, H2bA_controls, H2bB_controls,
                out= "Tables/H2control.tex",  style="qje",
                title="Testing H2a-H2b with Control variables",
                keep = "shaming",
                dep.var.labels = c("Favorability (Shamer)", "Military Coop. (Shamer)", 
                                   "Economic Coop. (Shamer)"),
                covariate.labels ="Shaming",
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:H2controls",
                # star.char = c("", "", ""),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{4}{c} {\\parbox[t]{14cm}{ \\textit{Notes:}
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex

tab1 <- gsub("\\begin{tabular}"
             ,"\\resizebox{1.1\\textwidth}{!}{\n\\begin{tabular}"
             ,tab1
             ,fixed=T)

tab1 <- gsub("\\end{tabular}"
             ,"\\end{tabular}}"
             ,tab1
             ,fixed=T)

cat(tab1, sep = "\n")

write(tab1, 'Tables/H2control.tex')



#Plot H1a, H1b, H2a, H2b
main_H12ab <- rbind(H1a, H1bA,H1bB,
                   H2a, H2bA,H2bB) %>% 
  filter(.,
         term == "shaming")

order_x <- c("Favorability \n (Target)", "Military Coop. \n (Target)", 
             "Economic Coop. \n (Target)", "Favorability \n (Shamer)", 
             "Military Coop. \n (Shamer)", "Economic Coop. \n (Shamer)") 

main_H12ab <- main_H12ab %>%
  mutate(color = ifelse(Towards == "shamer", "gray50", "black"))
colors <- main_H12ab$color[order(main_H12ab$Towards)]

# Generate main_H12ab plot 
ggplot(main_H12ab, aes(x = outcome, y = estimate)) +
  geom_hline(yintercept = 0, color = "gray50", linetype = 2, size = 0.2) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.6)) +
  geom_text(aes(label = round(estimate, digits = 2)),
            family = "Times", 
            position = position_dodge(width = 0.6),  # Align text with points
            hjust = -0.6,
            show_guide  = FALSE) + 
  scale_colour_grey()+
  scale_x_discrete(limits= order_x) +
  labs(x = "Outcome",
       y = "Effect Size of Shaming") +
  guides(text = "none") + 
  theme(text = element_text(size = 16, family = "Times"),
        legend.key=element_blank(),
        axis.text.x = element_text(color = c("gray50", "gray50", "gray50",
                                             "black","black", "black"), size = 14, margin = margin(0,0,20,0)),
        panel.grid.major = element_blank(),
        plot.caption = element_text(size = 14, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

ggsave("Figures/H12ab.pdf", width = 10, height = 5)


#Test H2c (Does the relationship with the shamer moderate H2a, H2bA, or H2bB?)

##favorability
H2c_Aally<-lm_robust(shamer_favorability~shaming,sham_ally)%>% tidy %>% 
  mutate(.,
         outcome = "Favorability \n (Shamer)",
         Towards = "shamer",
         Shamer = "U.S. Ally")

H2c_Aenemy<-lm_robust(shamer_favorability~shaming,sham_enemy)%>% tidy %>% 
  mutate(.,
         outcome = "Favorability \n (Shamer)",
         Towards = "shamer",
         Shamer = "U.S. Enemy")

##military
H2c_Bally<-lm_robust(shamer_military~shaming,sham_ally)%>% tidy %>% 
  mutate(.,
         outcome = "Military Coop. \n (Shamer)",
         Towards = "shamer",
         Shamer = "U.S. Ally")

H2c_Benemy<-lm_robust(shamer_military~shaming,sham_enemy)%>% tidy %>% 
  mutate(.,
         outcome = "Military Coop. \n (Shamer)",
         Towards = "shamer",
         Shamer = "U.S. Enemy")


##economic
H2c_Cally<-lm_robust(shamer_economic~shaming,sham_ally)%>% tidy %>% 
  mutate(.,
         outcome = "Economic Coop. \n (Shamer)",
         Towards = "shamer",
         Shamer = "U.S. Ally")

H2c_Cenemy<-lm_robust(shamer_economic~shaming,sham_enemy)%>% tidy %>% 
  mutate(.,
         outcome = "Economic Coop. \n (Shamer)",
         Towards = "shamer",
         Shamer = "U.S. Enemy")

#plot H1c and H2c results

#Plot H1a, H1b, H2a, H2b
main_H12c <- rbind(H1c_Aally, H1c_Aenemy,
                   H1c_Bally, H1c_Benemy,
                   H1c_Cally, H1c_Cenemy,
                   H2c_Aally, H2c_Aenemy,
                   H2c_Bally, H2c_Benemy,
                   H2c_Cally, H2c_Cenemy) %>% 
  filter(.,
         term == "shaming")

# Generate main_H12c plot 
ggplot(main_H12c, aes(x = outcome, y = estimate, colour = Shamer, shape = Shamer)) +
  geom_hline(yintercept = 0, color = "gray50", linetype = 2, size = 0.2) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.6)) +
  geom_text(aes(label = round(estimate, digits = 2)),
            family = "Times", 
            position = position_dodge(width = 0.6),  # Align text with points
            hjust = -0.6,
            show_guide  = FALSE) + 
  scale_color_manual(values = c("dodgerblue2","firebrick2")) +
  scale_fill_manual(values = c("dodgerblue2","firebrick2")) +
  #scale_colour_grey()+
  scale_x_discrete(limits= order_x) +
  labs(x = "Outcome",
       y = "Effect Size of Shaming") +
  guides(text = "none") + 
  theme(text = element_text(size = 16, family = "Times"),
        legend.key=element_blank(),
        axis.text.x = element_text(color = c("gray50", "gray50", "gray50",
                                             "black","black", "black"), size = 14, margin = margin(0,0,20,0)),
        panel.grid.major = element_blank(),
        plot.caption = element_text(size = 14, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

ggsave("Figures/H12c.pdf", width = 14, height = 5)

#report interaction terms of H2c in table form
##H2c_A:favorability
H2c_A<-lm(shamer_favorability~shaming*shamer_ally,data)

##H2c_B:military cooperation
H2c_B<-lm(shamer_military~shaming*shamer_ally,data)

##H2c_C:economic cooperation
H2c_C<-lm(shamer_economic~shaming*shamer_ally,data)


#with controls
H2c_Acontrols<-lm(shamer_favorability~shaming*shamer_ally+
                    Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state,data)


H2c_Bcontrols<-lm(shamer_military~shaming*shamer_ally+
                    Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state,data)


H2c_Ccontrols<-lm(shamer_economic~shaming*shamer_ally+
                    Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state,data)



#save table of interaction terms
#Create table
tab1<-stargazer(H1c_A, H1c_B, H1c_C,H2c_A, H2c_B, H2c_C,
                out= "Tables/H12c.tex",  style="qje",
                title="Heterogeneous treatmet effect of shaming*shamer ally",
                keep = c("shaming", "shamer_ally","shaming:shamer_ally"),
                dep.var.labels = c("Favorability (Target)", "Military Coop. (Target)", 
                                   "Economic Coop. (Target)",
                                   "Favorability (Shamer)","Military Coop. (Shamer)",
                                   "Economic Coop. (Shamer)"),
                covariate.labels = c("Shaming", "Shamer Ally", "Shaming*Shamer Ally"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:H12c",
                # star.char = c("", "", ""),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{7}{c} {\\parbox[t]{14cm}{ \\textit{Notes:}
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex

tab1 <- gsub("\\begin{tabular}"
               ,"\\resizebox{1.1\\textwidth}{!}{\n\\begin{tabular}"
               ,tab1
               ,fixed=T)

tab1 <- gsub("\\end{tabular}"
               ,"\\end{tabular}}"
               ,tab1
               ,fixed=T)

cat(tab1, sep = "\n")

write(tab1, 'Tables/H12c.tex')


#save table with controls
#Create table -- H1c controls
tab1<-stargazer(H1c_Acontrols, H1c_Bcontrols, 
                H1c_Ccontrols,
                out= "Tables/H1cControls.tex",  style="qje",
                title="Testing H1c with controls",
                keep = c("shaming", "shamer_ally","shaming:shamer_ally"),
                dep.var.labels = c("Favorability (Target)", "Military Coop. (Target)", 
                                   "Economic Coop. (Target)"),
                covariate.labels = c("Shaming", "Shamer Ally", "Shaming*Shamer Ally"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:H1cControls",
                # star.char = c("", "", ""),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{4}{c} {\\parbox[t]{14cm}{ \\textit{Notes:}
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex

tab1 <- gsub("\\begin{tabular}"
             ,"\\resizebox{1.1\\textwidth}{!}{\n\\begin{tabular}"
             ,tab1
             ,fixed=T)

tab1 <- gsub("\\end{tabular}"
             ,"\\end{tabular}}"
             ,tab1
             ,fixed=T)

cat(tab1, sep = "\n")

write(tab1, 'Tables/H1cControls.tex')


#Create table -- H2c controls
tab1<-stargazer(H2c_Acontrols, H2c_Bcontrols, H2c_Ccontrols,
                out= "Tables/H2cControls.tex",  style="qje",
                title="Testing H2c with controls",
                keep = c("shaming", "shamer_ally","shaming:shamer_ally"),
                dep.var.labels = c("Favorability (Shamer)","Military Coop. (Shamer)", 
                                   "Economic Coop. (Shamer)"),
                covariate.labels = c("Shaming", "Shamer Ally", "Shaming*Shamer Ally"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:H2cControls",
                # star.char = c("", "", ""),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{4}{c} {\\parbox[t]{14cm}{ \\textit{Notes:}
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex

tab1 <- gsub("\\begin{tabular}"
             ,"\\resizebox{1.1\\textwidth}{!}{\n\\begin{tabular}"
             ,tab1
             ,fixed=T)

tab1 <- gsub("\\end{tabular}"
             ,"\\end{tabular}}"
             ,tab1
             ,fixed=T)

cat(tab1, sep = "\n")

write(tab1, 'Tables/H2cControls.tex')



#Test H2d (Does the relationship with the shamer moderate H2a, H2bA, or 21bB?)

##favorability
H2d_Aally<-lm_robust(shamer_favorability~shaming,tar_ally)%>% tidy %>% 
  mutate(.,
         outcome = "Favorability \n (Shamer)",
         Towards = "shamer",
         Target = "U.S. Ally")

H2d_Aenemy<-lm_robust(shamer_favorability~shaming,tar_enemy)%>% tidy %>% 
  mutate(.,
         outcome = "Favorability \n (Shamer)",
         Towards = "shamer",
         Target = "U.S. Enemy") 


##military
H2d_Bally<-lm_robust(shamer_military~shaming,tar_ally)%>% tidy %>% 
  mutate(.,
         outcome = "Military Coop. \n (Shamer)",
         Towards = "shamer",
         Target = "U.S. Ally") 

H2d_Benemy<-lm_robust(shamer_military~shaming,tar_enemy)%>% tidy %>% 
  mutate(.,
         outcome = "Military Coop. \n (Shamer)",
         Towards = "shamer",
         Target = "U.S. Enemy") 


##economic
H2d_Cally<-lm_robust(shamer_economic~shaming,tar_ally)%>% tidy %>% 
  mutate(.,
         outcome = "Economic Coop. \n (Shamer)",
         Towards = "shamer",
         Target = "U.S. Ally") 

H2d_Cenemy<-lm_robust(shamer_economic~shaming,tar_enemy)%>% tidy %>% 
  mutate(.,
         outcome = "Economic Coop. \n (Shamer)",
         Towards = "shamer",
         Target = "U.S. Enemy") 



#plot H1c and H2c results

#Plot H1a, H1b, H2a, H2b
main_H12d <- rbind(H1d_Aally, H1d_Aenemy,
                   H1d_Bally, H1d_Benemy,
                   H1d_Cally, H1d_Cenemy,
                   H2d_Aally, H2d_Aenemy,
                   H2d_Bally, H2d_Benemy,
                   H2d_Cally, H2d_Cenemy) %>% 
  filter(.,
         term == "shaming")

# Generate main_H12d plot 
ggplot(main_H12d, aes(x = outcome, y = estimate, colour = Target, shape = Target)) +
  geom_hline(yintercept = 0, color = "gray50", linetype = 2, size = 0.2) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.6)) +
  geom_text(aes(label = round(estimate, digits = 2)),
            family = "Times", 
            position = position_dodge(width = 0.6),  # Align text with points
            hjust = -0.6,
            show_guide  = FALSE) + 
  scale_color_manual(values = c("dodgerblue2","firebrick2")) +
  scale_fill_manual(values = c("dodgerblue2","firebrick2")) +
  #scale_colour_grey()+
  scale_x_discrete(limits= order_x) +
  labs(x = "Outcome",
       y = "Effect Size of Shaming") +
  guides(text = "none") + 
  theme(text = element_text(size = 16, family = "Times"),
        legend.key=element_blank(),
        axis.text.x = element_text(color = c("gray50", "gray50", "gray50",
                                             "black","black", "black"), size = 14, margin = margin(0,0,20,0)),
        panel.grid.major = element_blank(),
        plot.caption = element_text(size = 14, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

ggsave("Figures/H12d.pdf", width = 14, height = 5)


#report interaction terms of H2d in table
##H2d_A:favorability
H2d_A<-lm(shamer_favorability~shaming*target_ally,data)

##H2d_B:military cooperation
H2d_B<-lm(shamer_military~shaming*target_ally,data)

##H2d_C:economic cooperation
H2d_C<-lm(shamer_economic~shaming*target_ally,data)

#with controls
H2d_Acontrols<-lm(shamer_favorability~shaming*target_ally+
                    Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state,data)

H2d_Bcontrols<-lm(shamer_military~shaming*target_ally+
                    Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state,data)

H2d_Ccontrols<-lm(shamer_economic~shaming*target_ally+
                    Democrat+Republican+Independent+
                    Educated+Black+White+Asian+Native+Other+Religious+
                    Female+Age+Political_interest+Conservative+hawk1+hawk2+
                    hawk3+hawk4+state,data)


#save table of interaction terms
#Create table
tab2<-stargazer(H1d_A, H1d_B, H1d_C,
                H2d_A, H2d_B, H2d_C,
                out= "Tables/H12d.tex",  style="qje",
                title="Heterogeneous treatmet effect of shaming*target ally",
                keep = c("shaming", "target_ally","shaming:target_ally"),
                dep.var.labels = c("Favorability (Target)", "Military Coop. (Target)", 
                                   "Economic Coop. (Target)",
                                   "Favorability (Shamer)","Military Coop. (Shamer)",
                                   "Economic Coop. (Shamer)"),
                covariate.labels = c("Shaming", "Target Ally", "Shaming*Target Ally"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:H12d",
                # star.char = c("", "", ""),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{7}{c} {\\parbox[t]{14cm}{ \\textit{Notes:}
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab2[grepl("Note",tab2)] <- note.latex

tab2 <- gsub("\\begin{tabular}"
             ,"\\resizebox{1.1\\textwidth}{!}{\n\\begin{tabular}"
             ,tab2
             ,fixed=T)

tab2 <- gsub("\\end{tabular}"
             ,"\\end{tabular}}"
             ,tab2
             ,fixed=T)

cat(tab2, sep = "\n")

write(tab2, 'Tables/H12d.tex')


#save table with controls
#Create table -- H1d controls
tab2<-stargazer(H1d_Acontrols, H1d_Bcontrols, H1d_Ccontrols,
                out= "Tables/H1dcontrols.tex",  style="qje",
                title="Testing H1d with Controls",
                keep = c("shaming", "target_ally","shaming:target_ally"),
                dep.var.labels = c("Favorability (Target)", "Military Coop. (Target)", 
                                   "Economic Coop. (Target)"),
                covariate.labels = c("Shaming", "Target Ally", "Shaming*Target Ally"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:H1dcontrols",
                # star.char = c("", "", ""),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{4}{c} {\\parbox[t]{14cm}{ \\textit{Notes:}
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab2[grepl("Note",tab2)] <- note.latex

tab2 <- gsub("\\begin{tabular}"
             ,"\\resizebox{1.1\\textwidth}{!}{\n\\begin{tabular}"
             ,tab2
             ,fixed=T)

tab2 <- gsub("\\end{tabular}"
             ,"\\end{tabular}}"
             ,tab2
             ,fixed=T)

cat(tab2, sep = "\n")

write(tab2, 'Tables/H1dcontrols.tex')


#save table with controls
#Create table -- H2d controls
tab2<-stargazer(H2d_Acontrols, H2d_Bcontrols, H2d_Ccontrols,
                out= "Tables/H2dcontrols.tex",  style="qje",
                title="Testing H2d with Controls",
                keep = c("shaming", "target_ally","shaming:target_ally"),
                dep.var.labels = c("Favorability (Shamer)","Military Coop. (Shamer)",
                                   "Economic Coop. (Shamer)"),
                covariate.labels = c("Shaming", "Target Ally", "Shaming*Target Ally"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:H2dcontrols",
                # star.char = c("", "", ""),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{4}{c} {\\parbox[t]{14cm}{ \\textit{Notes:}
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab2[grepl("Note",tab2)] <- note.latex

tab2 <- gsub("\\begin{tabular}"
             ,"\\resizebox{1.1\\textwidth}{!}{\n\\begin{tabular}"
             ,tab2
             ,fixed=T)

tab2 <- gsub("\\end{tabular}"
             ,"\\end{tabular}}"
             ,tab2
             ,fixed=T)

cat(tab2, sep = "\n")

write(tab2, 'Tables/H2dcontrols.tex')


#Test H3 (Does shaming decrease/increase support for international norms?)#

##H3a: detention
H3a<- lm_robust(detention_new~shaming, data = data)%>% tidy %>% 
  mutate(.,
         outcome = "Oppose \n Arbitrary Detention") 

##H3b: human rights
H3b<- lm_robust(HRs_new~shaming, data = data)%>% tidy %>% 
  mutate(.,
         outcome = "Support \n Human Rights") 

##H3c: human rights
H3c<- lm_robust(law_new~shaming, data = data)%>% tidy %>% 
  mutate(.,
         outcome = "Support \n International Law") 

##H3d: index
H3d<- lm_robust(H3_index~shaming, data = data)%>% tidy %>% 
  mutate(.,
         outcome = "Index") 

#add controls
H3acontrols<- lm(detention_new~shaming+
                   Democrat+Republican+Independent+
                   Educated+Black+White+Asian+Native+Other+Religious+
                   Female+Age+Political_interest+Conservative+hawk1+hawk2+
                   hawk3+hawk4+state, data) 

H3bcontrols<- lm(HRs_new~shaming+
                   Democrat+Republican+Independent+
                   Educated+Black+White+Asian+Native+Other+Religious+
                   Female+Age+Political_interest+Conservative+hawk1+hawk2+
                   hawk3+hawk4+state, data) 

H3ccontrols<- lm(law_new~shaming+
                   Democrat+Republican+Independent+
                   Educated+Black+White+Asian+Native+Other+Religious+
                   Female+Age+Political_interest+Conservative+hawk1+hawk2+
                   hawk3+hawk4+state, data) 

H3dcontrols<- lm(H3_index~shaming+
                   Democrat+Republican+Independent+
                   Educated+Black+White+Asian+Native+Other+Religious+
                   Female+Age+Political_interest+Conservative+hawk1+hawk2+
                   hawk3+hawk4+state, data) 

summary(H3bcontrols)

#Plot H3
main_H3 <- rbind(H3a, H3b, H3c, H3d) %>% 
  filter(.,
         term == "shaming")

order_x <- c("Oppose \n Arbitrary Detention", "Support \n Human Rights", 
             "Support \n International Law", "Index") 

# Generate coefficient plot 
ggplot(main_H3, aes(x = outcome, y = estimate)) +
  geom_hline(yintercept = 0, color = "gray50", linetype = 2, size = 0.2) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.4),
                  fill = "white") +
  geom_text(aes(label= round(estimate, digits = 2)), hjust=-.7, size = 3,
            family = "Times") +
  #geom_text(aes(label= paste("(",round(std.error, digits = 1),")")), hjust=-.2, vjust =2 ,  size = 3, 
    #        family = "Times") +
  #coord_flip()+
  labs(x = "",
       y = "Effect Size of Shaming") +
  scale_x_discrete(limits= order_x) +
  theme(text = element_text(size = 10, family = "Times"),
        legend.key=element_blank(),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(size = 10),
        plot.caption = element_text(size = 10, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

ggsave("Figures/H3.pdf", width = 10, height = 5)

#save table with controls

#save table with controls
#Create table -- H3 controls
tab2<-stargazer(H3acontrols, H3bcontrols, H3ccontrols, H3dcontrols,
                out= "Tables/H3controls.tex",  style="qje",
                title="Testing H3 with Controls",
                keep = "shaming",
                dep.var.labels = c("Oppose Arbitrary Detention", "Support HRs", 
                                   "Support Int. Law", "Index"),
                covariate.labels = "Shaming",
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:H3controls",
                # star.char = c("", "", ""),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{5}{c} {\\parbox[t]{14cm}{ \\textit{Notes:}
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab2[grepl("Note",tab2)] <- note.latex

tab2 <- gsub("\\begin{tabular}"
             ,"\\resizebox{1.1\\textwidth}{!}{\n\\begin{tabular}"
             ,tab2
             ,fixed=T)

tab2 <- gsub("\\end{tabular}"
             ,"\\end{tabular}}"
             ,tab2
             ,fixed=T)

cat(tab2, sep = "\n")

write(tab2, 'Tables/H3controls.tex')


## ----------------------------------------------------------------------
# In-text discussion
# ----------------------------------------------------------------------

not_shamed <- data[ which(data$shaming==0), ]#control respondents (baseline)

#mean attitude towards ally-shamers, control
not_shamed_ally <- not_shamed[ which(not_shamed$shamer_ally==1), ] 
mean(not_shamed_ally$shamer_favorability,na.rm=T)

#mean attitude towards enemy-targets, control
not_shamed_enemy <- not_shamed[ which(not_shamed$target_ally==0), ]
mean(not_shamed_enemy$target_favorability,na.rm=T)



## ----------------------------------------------------------------------
# Supplementary Analyses
# ----------------------------------------------------------------------

#Decriptive Statistics


disc_table <- dplyr::select(data, c(target_favorability, shamer_favorability, 
                                    target_military, shamer_military, 
                                    target_economic, shamer_economic,
                                    detention_new, HRs_new, law_new, H3_index,
                                    Democrat, Republican,Independent, Educated, 
                                    Black,White, Asian, Native, Other, Religious, Female, 
                                    Political_interest, Hawkishness, Conservative, Age))

#save to Tables
stargazer(as.data.frame(disc_table),
          covariate.labels = c("Target Favorability", "Shamer Favorability", 
                               "Target military coop", "Shamer military coop", 
                               "Target economic coop", "Shamer economic coop",
                               "Oppose Detention", "Support HRs", "Support Int. law", "H3 Index",
                               "Democrat","Republican","Independent", "Educated", 
                               "Black","White", "Asian", "Native", "Other", "Religious", "Female", 
                               "Political interest", "Hawkishness", "Conservative", "Age"),
          
          title = "Descriptive Statistics",
          label = "tab:desc",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          type = "latex", out = "Tables/descriptives.tex")


#How are individuals classifying the relationship (ally/enemy) between US and the countries?
data_complete<-data[ which(!is.na(data$state)), ]
df_long <- data_complete[c("man3_1","man3_2","man3_3","man3_4","man3_5","man3_6","man3_7","man3_8","man3_9","man3_10")]
df_long <- df_long %>%
  pivot_longer(cols = everything(), names_to = "Country", values_to = "Response") %>%
  mutate(Response = factor(Response, levels = c("Enemy", "Unfriendly", "Friendly", "Ally"), labels = c(1, 2, 3, 4)),
         Country = case_when(
           Country=="man3_1"~"Australia", Country=="man3_2"~"Canada", Country=="man3_3"~"France", Country=="man3_4"~"Japan",
           Country=="man3_5"~"Israel", Country=="man3_6"~"China",Country=="man3_7"~"Iran",Country=="man3_8"~"North Korea",
           Country=="man3_9"~"Russia",Country=="man3_10"~"Venezuela"
         ))

df_avg <- df_long %>%
  group_by(Country) %>%
  summarize(AverageResponse = mean(as.numeric(Response),na.rm=T))

# Plot the average response for each country
ggplot(df_avg, aes(x = reorder(Country, AverageResponse), y = AverageResponse, fill = AverageResponse)) +
  geom_bar(stat = "identity", width = 0.7) + 
  scale_fill_gradient(low = "#D3D3D3", high = "#808080") +  
  labs(x = "Country", y = "Average score", 
       title = "Perceptions of U.S. relationship with each country",
       subtitle = "Scores range from Enemy (1) to Ally (4)") +
  theme_minimal(base_size = 10) + 
  geom_text(aes(label = round(AverageResponse, 1)), vjust = -0.5, size = 3, family = "Times") +
# Set base font size
  theme(
    text = element_text(size = 10, family = "Times"),
    legend.key=element_blank(),
    panel.grid.major = element_blank(), 
    axis.text.x = element_text(size = 10),
    plot.caption = element_text(size = 10, family = "Times",hjust = -.02),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(), 
    legend.position = "none" ,
    axis.line = element_line(colour = "black")) 

ggsave("Figures/ally_perception.pdf", width = 10, height = 5)

#Do perceptions of shamer's ally status mediate the effect of shaming? (H1c, H2c)

##H1c mediated -- attitudes toward *target*

H1c_med <- lm(shamer_allyperception~shaming*shamer_ally,data=data)

#Outcome 1 (favorability)
H1cA_out <- lm(target_favorability~shaming*shamer_ally + shamer_allyperception*shamer_ally, data = data)
summary(H1cA_out)
# Test the mediation effect at different levels of ally status (ally = 0 and ally = 1)
mediation1_ally0 <- mediate(H1c_med, H1cA_out, treat = "shaming", 
                                mediator = "shamer_allyperception", covariates = list(shamer_ally = 0), sims = 100)

mediation1_ally1 <- mediate(H1c_med, H1cA_out, treat = "shaming", 
                            mediator = "shamer_allyperception", covariates = list(shamer_ally = 1), sims = 100)


#plot and save png
png("Figures/mediation1.png")
par(mfrow = c(1, 2))

plot(mediation1_ally0, main = "Shamer = Enemy", labels = c("ACME\n(Ally \n perception)",
                                                                                "Direct \nEffect \n(Target \n favorability)", "Total \nEffect"),
     col = "blue", cex.axis = 0.6)

plot(mediation1_ally1, main = "Shamer = Ally", labels = c("ACME\n(Ally \n perception)",
                                                                                "Direct \nEffect \n(Target \n favorability)", "Total \nEffect"),
     col = "red", cex.axis = 0.6)

dev.off()


#Outcome 2 (military)
H1cB_out <- lm(target_military~shaming*shamer_ally + shamer_allyperception*shamer_ally, data = data)

# Test the mediation effect at different levels of ally status (ally = 0 and ally = 1)
mediation2_ally0 <- mediate(H1c_med, H1cB_out, treat = "shaming", 
                            mediator = "shamer_allyperception", covariates = list(shamer_ally = 0), sims = 100)

mediation2_ally1 <- mediate(H1c_med, H1cB_out, treat = "shaming", 
                            mediator = "shamer_allyperception", covariates = list(shamer_ally = 1), sims = 100)


#plot and save png
png("Figures/mediation2.png")
par(mfrow = c(1, 2))

plot(mediation2_ally0, main = "Shamer = Enemy", labels = c("ACME\n(Ally \n perception)",
                                                           "Direct \nEffect \n(Target \n military)", "Total \nEffect"),
     col = "blue", cex.axis = 0.6)

plot(mediation2_ally1, main = "Shamer = Ally", labels = c("ACME\n(Ally \n perception)",
                                                          "Direct \nEffect \n(Target \n military)", "Total \nEffect"),
     col = "red", cex.axis = 0.6)

dev.off()



#Outcome 3 (economic)
H1cC_out <- lm(target_economic~shaming*shamer_ally + shamer_allyperception*shamer_ally, data = data)

# Test the mediation effect at different levels of ally status (ally = 0 and ally = 1)
mediation3_ally0 <- mediate(H1c_med, H1cC_out, treat = "shaming", 
                            mediator = "shamer_allyperception", covariates = list(shamer_ally = 0), sims = 100)

mediation3_ally1 <- mediate(H1c_med, H1cC_out, treat = "shaming", 
                            mediator = "shamer_allyperception", covariates = list(shamer_ally = 1), sims = 100)


#plot and save png
png("Figures/mediation3.png")
par(mfrow = c(1, 2))

plot(mediation3_ally0, main = "Shamer = Enemy", labels = c("ACME\n(Ally \n perception)",
                                                           "Direct \nEffect \n(Target \n economic)", "Total \nEffect"),
     col = "blue", cex.axis = 0.6)

plot(mediation3_ally1, main = "Shamer = Ally", labels = c("ACME\n(Ally \n perception)",
                                                          "Direct \nEffect \n(Target \n economic)", "Total \nEffect"),
     col = "red", cex.axis = 0.6)

dev.off()


##H2c mediated -- attitudes toward *shamer*
H2c_med <- lm(shamer_allyperception~shaming*shamer_ally,data=data)

#Outcome 1 (favorability)
H2cA_out <- lm(shamer_favorability~shaming*shamer_ally + shamer_allyperception*shamer_ally, data = data)

# Test the mediation effect at different levels of ally status (ally = 0 and ally = 1)
mediation4_ally0 <- mediate(H2c_med, H2cA_out, treat = "shaming", 
                            mediator = "shamer_allyperception", covariates = list(shamer_ally = 0), sims = 100)

mediation4_ally1 <- mediate(H2c_med, H2cA_out, treat = "shaming", 
                            mediator = "shamer_allyperception", covariates = list(shamer_ally = 1), sims = 100)


#plot and save png
png("Figures/mediation4.png")
par(mfrow = c(1, 2))

plot(mediation4_ally0, main = "Shamer = Enemy", labels = c("ACME\n(Ally \n perception)",
                                                           "Direct \nEffect \n(Shamer \n favorability)", "Total \nEffect"),
     col = "blue", cex.axis = 0.6)

plot(mediation4_ally1, main = "Shamer = Ally", labels = c("ACME\n(Ally \n perception)",
                                                          "Direct \nEffect \n(Shamer \n favorability)", "Total \nEffect"),
     col = "red", cex.axis = 0.6)

dev.off()


#Outcome 2 (military)
H2cB_out <- lm(shamer_military~shaming*shamer_ally + shamer_allyperception*shamer_ally, data = data)

# Test the mediation effect at different levels of ally status (ally = 0 and ally = 1)
mediation5_ally0 <- mediate(H2c_med, H2cB_out, treat = "shaming", 
                            mediator = "shamer_allyperception", covariates = list(shamer_ally = 0), sims = 100)

mediation5_ally1 <- mediate(H2c_med, H2cB_out, treat = "shaming", 
                            mediator = "shamer_allyperception", covariates = list(shamer_ally = 1), sims = 100)


#plot and save png
png("Figures/mediation5.png")
par(mfrow = c(1, 2))

plot(mediation5_ally0, main = "Shamer = Enemy", labels = c("ACME\n(Ally \n perception)",
                                                           "Direct \nEffect \n(Shamer \n military)", "Total \nEffect"),
     col = "blue", cex.axis = 0.6)

plot(mediation5_ally1, main = "Shamer = Ally", labels = c("ACME\n(Ally \n perception)",
                                                          "Direct \nEffect \n(Shamer \n military)", "Total \nEffect"),
     col = "red", cex.axis = 0.6)

dev.off()



#Outcome 3 (economic)
H2cC_out <- lm(shamer_economic~shaming*shamer_ally + shamer_allyperception*shamer_ally, data = data)

# Test the mediation effect at different levels of ally status (ally = 0 and ally = 1)
mediation6_ally0 <- mediate(H2c_med, H2cC_out, treat = "shaming", 
                            mediator = "shamer_allyperception", covariates = list(shamer_ally = 0), sims = 100)

mediation6_ally1 <- mediate(H2c_med, H2cC_out, treat = "shaming", 
                            mediator = "shamer_allyperception", covariates = list(shamer_ally = 1), sims = 100)


#plot and save png
png("Figures/mediation6.png")
par(mfrow = c(1, 2))

plot(mediation6_ally0, main = "Shamer = Enemy", labels = c("ACME\n(Ally \n perception)",
                                                           "Direct \nEffect \n(Shamer \n economic)", "Total \nEffect"),
     col = "blue", cex.axis = 0.6)

plot(mediation6_ally1, main = "Shamer = Ally", labels = c("ACME\n(Ally \n perception)",
                                                          "Direct \nEffect \n(Shamer \n economic)", "Total \nEffect"),
     col = "red", cex.axis = 0.6)

dev.off()


#Do perceptions of target's ally status mediate the effect of shaming? (H1c, H2c)

##H1d mediated -- attitudes toward *target* (target ally)

H1d_med <- lm(target_allyperception~shaming*target_ally,data=data)

#Outcome 1 (favorability)
H1dA_out <- lm(target_favorability~shaming*target_ally + target_allyperception*target_ally, data = data)
summary(H1dA_out)
# Test the mediation effect at different levels of ally status (ally = 0 and ally = 1)
mediation7_ally0 <- mediate(H1d_med, H1dA_out, treat = "shaming", 
                            mediator = "target_allyperception", covariates = list(target_ally = 0), sims = 100)

mediation7_ally1 <- mediate(H1d_med, H1dA_out, treat = "shaming", 
                            mediator = "target_allyperception", covariates = list(target_ally = 1), sims = 100)


#plot and save png
png("Figures/mediation7.png")
par(mfrow = c(1, 2))

plot(mediation7_ally0, main = "Target = Enemy", labels = c("ACME\n(Ally \n perception)",
                                                           "Direct \nEffect \n(Target \n favorability)", "Total \nEffect"),
     col = "blue", cex.axis = 0.6)

plot(mediation7_ally1, main = "Target = Ally", labels = c("ACME\n(Ally \n perception)",
                                                          "Direct \nEffect \n(Target \n favorability)", "Total \nEffect"),
     col = "red", cex.axis = 0.6)

dev.off()


#Outcome 2 (military)
H1dB_out <- lm(target_military~shaming*target_ally + target_allyperception*target_ally, data = data)

# Test the mediation effect at different levels of ally status (ally = 0 and ally = 1)
mediation8_ally0 <- mediate(H1d_med, H1dB_out, treat = "shaming", 
                            mediator = "target_allyperception", covariates = list(target_ally = 0), sims = 100)

mediation8_ally1 <- mediate(H1d_med, H1dB_out, treat = "shaming", 
                            mediator = "target_allyperception", covariates = list(target_ally = 1), sims = 100)


#plot and save png
png("Figures/mediation8.png")
par(mfrow = c(1, 2))

plot(mediation8_ally0, main = "Target = Enemy", labels = c("ACME\n(Ally \n perception)",
                                                           "Direct \nEffect \n(Target \n military)", "Total \nEffect"),
     col = "blue", cex.axis = 0.6)

plot(mediation8_ally1, main = "Target = Ally", labels = c("ACME\n(Ally \n perception)",
                                                          "Direct \nEffect \n(Target \n military)", "Total \nEffect"),
     col = "red", cex.axis = 0.6)

dev.off()



#Outcome 3 (economic)
H1dC_out <- lm(target_economic~shaming*target_ally + target_allyperception*target_ally, data = data)

# Test the mediation effect at different levels of ally status (ally = 0 and ally = 1)
mediation9_ally0 <- mediate(H1d_med, H1dC_out, treat = "shaming", 
                            mediator = "target_allyperception", covariates = list(target_ally = 0), sims = 100)

mediation9_ally1 <- mediate(H1d_med, H1dC_out, treat = "shaming", 
                            mediator = "target_allyperception", covariates = list(target_ally = 1), sims = 100)


#plot and save png
png("Figures/mediation9.png")
par(mfrow = c(1, 2))

plot(mediation9_ally0, main = "Target = Enemy", labels = c("ACME\n(Ally \n perception)",
                                                           "Direct \nEffect \n(Target \n economic)", "Total \nEffect"),
     col = "blue", cex.axis = 0.6)

plot(mediation9_ally1, main = "Target = Ally", labels = c("ACME\n(Ally \n perception)",
                                                          "Direct \nEffect \n(Target \n economic)", "Total \nEffect"),
     col = "red", cex.axis = 0.6)

dev.off()


##H2d mediated -- attitudes toward *shamer*
H2d_med <- lm(target_allyperception~shaming*target_ally,data=data)

#Outcome 1 (favorability)
H2dA_out <- lm(shamer_favorability~shaming*target_ally + target_allyperception*target_ally, data = data)

# Test the mediation effect at different levels of ally status (ally = 0 and ally = 1)
mediation10_ally0 <- mediate(H2d_med, H2dA_out, treat = "shaming", 
                            mediator = "target_allyperception", covariates = list(target_ally = 0), sims = 100)

mediation10_ally1 <- mediate(H2d_med, H2dA_out, treat = "shaming", 
                            mediator = "target_allyperception", covariates = list(target_ally = 1), sims = 100)


#plot and save png
png("Figures/mediation10.png")
par(mfrow = c(1, 2))

plot(mediation10_ally0, main = "Target = Enemy", labels = c("ACME\n(Ally \n perception)",
                                                           "Direct \nEffect \n(Shamer \n favorability)", "Total \nEffect"),
     col = "blue", cex.axis = 0.6)

plot(mediation10_ally1, main = "Target = Ally", labels = c("ACME\n(Ally \n perception)",
                                                          "Direct \nEffect \n(Shamer \n favorability)", "Total \nEffect"),
     col = "red", cex.axis = 0.6)

dev.off()


#Outcome 2 (military)
H2dB_out <- lm(shamer_military~shaming*target_ally + target_allyperception*target_ally, data = data)

# Test the mediation effect at different levels of ally status (ally = 0 and ally = 1)
mediation11_ally0 <- mediate(H2d_med, H2dB_out, treat = "shaming", 
                            mediator = "target_allyperception", covariates = list(target_ally = 0), sims = 100)

mediation11_ally1 <- mediate(H2d_med, H2dB_out, treat = "shaming", 
                            mediator = "target_allyperception", covariates = list(target_ally = 1), sims = 100)


#plot and save png
png("Figures/mediation11.png")
par(mfrow = c(1, 2))

plot(mediation11_ally0, main = "Target = Enemy", labels = c("ACME\n(Ally \n perception)",
                                                           "Direct \nEffect \n(Shamer \n military)", "Total \nEffect"),
     col = "blue", cex.axis = 0.6)

plot(mediation11_ally1, main = "Target = Ally", labels = c("ACME\n(Ally \n perception)",
                                                          "Direct \nEffect \n(Shamer \n military)", "Total \nEffect"),
     col = "red", cex.axis = 0.6)

dev.off()



#Outcome 3 (economic)
H2dC_out <- lm(shamer_economic~shaming*target_ally + target_allyperception*target_ally, data = data)

# Test the mediation effect at different levels of ally status (ally = 0 and ally = 1)
mediation12_ally0 <- mediate(H2d_med, H2dC_out, treat = "shaming", 
                            mediator = "target_allyperception", covariates = list(target_ally = 0), sims = 100)

mediation12_ally1 <- mediate(H2d_med, H2dC_out, treat = "shaming", 
                            mediator = "target_allyperception", covariates = list(target_ally = 1), sims = 100)


#plot and save png
png("Figures/mediation12.png")
par(mfrow = c(1, 2))

plot(mediation12_ally0, main = "Target = Enemy", labels = c("ACME\n(Ally \n perception)",
                                                           "Direct \nEffect \n(Shamer \n economic)", "Total \nEffect"),
     col = "blue", cex.axis = 0.6)

plot(mediation12_ally1, main = "Target = Ally", labels = c("ACME\n(Ally \n perception)",
                                                          "Direct \nEffect \n(Shamer \n economic)", "Total \nEffect"),
     col = "red", cex.axis = 0.6)

dev.off()




#How are individuals classifying the relationship (allies/enemies) between country-pairs?

# Calculate the average score for each shamer-target pair
df_avg_scores <- data_complete %>%
  group_by(shamer_name, target_name) %>%
  summarise(avg_score = mean(man4_numeric, na.rm = TRUE)) %>%
  ungroup()

order_xy<-c("Australia" , "Canada"    ,   "France"    ,     
           "Israel"   ,   "Japan"   , "China" ,  "Iran" ,   "North Korea" , "Russia"   ,   "Venezuela"  )


# Create a heatmap of the average scores
ggplot(df_avg_scores, aes(x = target_name, y = shamer_name, fill = avg_score)) +
  geom_tile() +
  scale_fill_gradient(low = "lightgray", high = "black") +
  labs(x = "Target", y = "Shamer", fill = "Average Score",
       title = "Perceptions of relationship between targets and shamers",
       subtitle = "Scores range from Enemies (1) to Allies (4)") +
  theme_minimal() +
  scale_x_discrete(limits= order_xy) +
  scale_y_discrete(limits= order_xy) +
  theme(
        text = element_text(size = 10, family = "Times"),
        legend.key=element_blank(),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(angle = 45, hjust = 1),
        plot.caption = element_text(size = 10, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black")) 

ggsave("Figures/ally_perception2.pdf", width = 6, height = 5)


#Do perceptions of shamer-target ally status moderate the effect of shaming?
##When pair is classified as allies, is shaming more likely to increase support for shamer?
##Sup1_A:favorability
Sup1_A<-lm(shamer_favorability~shaming*man4_numeric,data)
Sup1_a<-lm(shamer_favorability~shaming*man4_numeric+shaming*shamer_ally,data)

##Sup1_B:military cooperation
Sup1_B<-lm(shamer_military~shaming*man4_numeric,data)
Sup1_b<-lm(shamer_military~shaming*man4_numeric+shaming*shamer_ally,data)

##Sup1_C:economic cooperation
Sup1_C<-lm(shamer_economic~shaming*man4_numeric,data)
Sup1_c<-lm(shamer_economic~shaming*man4_numeric+shaming*shamer_ally,data)

#save table of interaction terms
#Create table
tab5<-stargazer(Sup1_A,Sup1_a,
                Sup1_B, Sup1_b, 
                Sup1_C, Sup1_c,
                out= "Tables/Sup1.tex",  style="qje",
                title="Moderating effect of Target-Shamer Ally Status (1)",
                keep = c("shaming", "man4_numeric","shamer_ally", "shaming:man4_numeric",
                         "shaming:shamer_ally"),
                dep.var.labels = c("Shamer Favorability", "Military Cooperation w. Shamer", "Economic Cooperation w. Shamer"),
                covariate.labels = c("Shaming", "Allies", "Shamer ally", "Shaming*Allies",
                                     "Shaming*Ally"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = c(0.05, 0.01, 0.001),
                label = "tab:sup1",
                #star.char = c("", "", "", ""),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{7}{c} {\\parbox[t]{12cm}{ \\textit{Notes:}
* $p < .05$, ** $p < .01$, *** $p < .001$}} \\\\"
tab5[grepl("Note",tab5)] <- note.latex
tab5 <- gsub("\\begin{tabular}"
             ,"\\resizebox{1.1\\textwidth}{!}{\n\\begin{tabular}"
             ,tab5
             ,fixed=T)

tab5 <- gsub("\\end{tabular}"
             ,"\\end{tabular}}"
             ,tab5
             ,fixed=T)

cat(tab5, sep = "\n")
write(tab5, 'Tables/Sup1.tex')

##When pair is classified as allies, is shaming more likely to decrease support for target?
##Sup2_A:favorability
Sup2_A<-lm(target_favorability~shaming*man4_numeric,data)

##Sup2_B:military cooperation
Sup2_B<-lm(target_military~shaming*man4_numeric,data)

##Sup2_C:economic cooperation
Sup2_C<-lm(target_economic~shaming*man4_numeric,data)

#save table of interaction terms
#Create table
tab6<-stargazer(Sup2_A, Sup2_B, Sup2_C,
                out= "Tables/Sup2.tex",  style="qje",
                title="Moderating effect of Target-Shamer Ally Status (2)",
                keep = c("shaming", "man4_numeric","shaming:man4_numeric"),
                dep.var.labels = c("Target Favorability", "Military Cooperation w. Target", "Economic Cooperation w. Target"),
                covariate.labels = c("Shaming", "Allies", "Shaming*Allies"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = c(0.05, 0.01, 0.001),
                label = "tab:sup2",
                # star.char = c("", "", ""),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))


note.latex <- "\\multicolumn{4}{c} {\\parbox[t]{12cm}{ \\textit{Notes:}
* $p < .05$, ** $p < .01$, *** $p < .001$}} \\\\"
tab6[grepl("Note",tab6)] <- note.latex
tab6 <- gsub("\\begin{tabular}"
             ,"\\resizebox{1.1\\textwidth}{!}{\n\\begin{tabular}"
             ,tab6
             ,fixed=T)

tab6 <- gsub("\\end{tabular}"
             ,"\\end{tabular}}"
             ,tab6
             ,fixed=T)

cat(tab6, sep = "\n")
write(tab6, 'Tables/Sup2.tex')

#HTEs -- Does treatment effect vary by respondents’ partisan affiliation (Democrat, Republican, Independent)?
#democrats
dem1<-data%>%filter(Democrat==1)%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Model = "Democrats") 

dem2<-data%>%filter(Democrat==1)%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Model = "Democrats") 

dem3<-data%>%filter(Democrat==1)%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Model = "Democrats") 

dem4<-data%>%filter(Democrat==1)%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Model = "Democrats") 

dem5<-data%>%filter(Democrat==1)%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Model = "Democrats") 

dem6<-data%>%filter(Democrat==1)%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Model = "Democrats") 

dem7<-data%>%filter(Democrat==1)%>%
  lm_robust(detention_new~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Oppose Detention",
         Model = "Democrats") 

dem8<-data%>%filter(Democrat==1)%>%
  lm_robust(HRs_new~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Support HRs",
         Model = "Democrats") 

dem9<- data%>%filter(Democrat==1)%>%
  lm_robust(law_new~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Support Int Law",
         Model = "Democrats") 

#republicans
rep1<-data%>%filter(Republican==1)%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Model = "Republicans") 

rep2<-data%>%filter(Republican==1)%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Model = "Republicans") 

rep3<-data%>%filter(Republican==1)%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Model = "Republicans") 

rep4<-data%>%filter(Republican==1)%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Model = "Republicans") 

rep5<-data%>%filter(Republican==1)%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Model = "Republicans") 

rep6<-data%>%filter(Republican==1)%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Model = "Republicans") 

rep7<-data%>%filter(Republican==1)%>%
  lm_robust(detention_new~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Oppose Detention",
         Model = "Republicans") 

rep8<-data%>%filter(Republican==1)%>%
  lm_robust(HRs_new~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Support HRs",
         Model = "Republicans") 

rep9<- data%>%filter(Republican==1)%>%
  lm_robust(law_new~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Support Int Law",
         Model = "Republicans") 

#independents
ind1<-data%>%filter(Independent==1)%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Model = "Independents") 

ind2<-data%>%filter(Independent==1)%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Model = "Independents") 

ind3<-data%>%filter(Independent==1)%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Model = "Independents") 

ind4<-data%>%filter(Independent==1)%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Model = "Independents") 

ind5<-data%>%filter(Independent==1)%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Model = "Independents") 

ind6<-data%>%filter(Independent==1)%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Model = "Independents") 

ind7<-data%>%filter(Independent==1)%>%
  lm_robust(detention_new~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Oppose Detention",
         Model = "Independents") 

ind8<-data%>%filter(Independent==1)%>%
  lm_robust(HRs_new~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Support HRs",
         Model = "Independents") 

ind9<- data%>%filter(Independent==1)%>%
  lm_robust(law_new~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Support Int Law",
         Model = "Independents") 

order_x <- c("Democrats", "Independents", "Republicans") 

#plot
model_party <- rbind(dem1,dem2,dem3,dem4,dem5,dem6,dem7,dem8,dem9,
                     ind1,ind2,ind3,ind4,ind5,ind6,ind7,ind8,ind9,
                     rep1,rep2,rep3,rep4,rep5,rep6,rep7,rep8,rep9) %>% 
  filter(.,
         term == "shaming")

model_party$Outcome <- factor(model_party$Outcome, levels = c("Target Favorability", 
                                                                            "Military Cooperation w. Target", 
                                                                            "Economic Cooperation w. Target",
                                                                            "Shamer Favorability",
                                                                            "Military Cooperation w. Shamer", 
                                                                            "Economic Cooperation w. Shamer",
                                                              "Oppose Detention",
                                                              "Support HRs","Support Int Law"))

# Generate coefficient plot 
ggplot(model_party, aes(x = Outcome, y = estimate, color = Model, shape = Model)) +
  geom_hline(yintercept = 0, color = "gray50", linetype = 2, size = 0.2) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.4)) +
  scale_colour_grey()+
  #scale_shape_manual(values=seq(0,9))+
  labs(x = "",
       y = "Effect Size of Shaming") +
  theme(text = element_text(size = 16, family = "Times"),
        legend.key=element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 16),
        panel.grid.major = element_blank(),
        plot.caption = element_text(size = 13, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

ggsave("Figures/byparty.pdf", width = 15, height = 5)


#Manipulation check:report proportion of respondents who passed shaming manipulation check

# Compute proportions for the different groups and conditions
shamed_respondents <- data[ which(data$shaming==1), ]
overall_man1_prop <- sum(data$man1_correct, na.rm=T) / nrow(na.omit(as.data.frame(data$man1)))
overall_man2_prop <- sum(data$man2_correct, na.rm=T) / nrow(na.omit(as.data.frame(data$man2)))
nonshamed_respondents <- data[ which(data$shaming==0), ]
shamed_man1_prop <- sum(shamed_respondents$man1_correct, na.rm=T) / nrow(na.omit(as.data.frame(shamed_respondents$man1)))
nonshamed_man1_prop <- sum(nonshamed_respondents$man1_correct, na.rm=T) / nrow(na.omit(as.data.frame(nonshamed_respondents$man1)))

shamed_man2_prop <- sum(shamed_respondents$man2_correct, na.rm=T) / nrow(na.omit(as.data.frame(shamed_respondents$man2)))
nonshamed_man2_prop <- sum(nonshamed_respondents$man2_correct, na.rm=T) / nrow(na.omit(as.data.frame(nonshamed_respondents$man2)))

# Create a data frame for plotting
manipulation_check_data <- data.frame(
  Condition = rep(c("Overall", "Shamed", "Non-shamed"), each = 2),
  Manipulation_Check = rep(c("Man1", "Man2"), times = 3),
  Proportion = c(overall_man1_prop, overall_man2_prop,
                 shamed_man1_prop, shamed_man2_prop,
                 nonshamed_man1_prop, nonshamed_man2_prop)
)

manipulation_check_data_filtered <- manipulation_check_data[ 
  !(manipulation_check_data$Manipulation_Check == "Man2" & manipulation_check_data$Condition != "Shamed"), 
]

order_x <- c("Overall", "Non-shamed", "Shamed") 

ggplot(manipulation_check_data_filtered, aes(x = Condition, y = Proportion, fill = Manipulation_Check)) +
  geom_bar(stat = "identity", position = "dodge") +
  theme_minimal() +
  labs(
    y = "Proportion Passed", 
    title = "Manipulation Check Pass Rates by Condition"
  ) +
  scale_fill_manual(
    values = c("darkgray", "lightgray"),
    labels = c("Recalled shaming", "Recalled shamer")
  ) +
  guides(fill = guide_legend(title = "Manipulation Checks")) + 
  ylim(0,1)+
  scale_x_discrete(limits= order_x) +
  theme(text = element_text(size = 10, family = "Times"),
        panel.grid.major = element_blank(),
        plot.caption = element_text(size = 10, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black")
  )

ggsave("Figures/manipulation.pdf", width = 6, height = 6)



#Ceiling/floor effects
#Drop one country at a time to examine whether one country is driving the effect

#Drop one *target* at a time to see if driving effect -- 
#Drop Australia
australia_tar1<-data%>%filter(target_name!="Canada")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Target = "Drop Australia") 

australia_tar2<-data%>%filter(target_name!="Australia")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Target = "Drop Australia") 

australia_tar3<-data%>%filter(target_name!="Australia")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Target = "Drop Australia") 

australia_tar4<-data%>%filter(target_name!="Australia")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Target = "Drop Australia") 

australia_tar5<-data%>%filter(target_name!="Australia")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Target = "Drop Australia") 

australia_tar6<-data%>%filter(target_name!="Australia")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Target = "Drop Australia") 

#Canada
canada_tar1<-data%>%filter(target_name!="Canada")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Target = "Drop Canada") 

canada_tar2<-data%>%filter(target_name!="Canada")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Target = "Drop Canada") 

canada_tar3<-data%>%filter(target_name!="Canada")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Target = "Drop Canada") 


canada_tar4<-data%>%filter(target_name!="Canada")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Target = "Drop Canada") 

canada_tar5<-data%>%filter(target_name!="Canada")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Target = "Drop Canada") 

canada_tar6<-data%>%filter(target_name!="Canada")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Target = "Drop Canada") 


#France
fance_tar1<-data%>%filter(target_name!="France")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Target = "Drop France") 

fance_tar2<-data%>%filter(target_name!="France")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Target = "Drop France") 

fance_tar3<-data%>%filter(target_name!="France")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Target = "Drop France") 

fance_tar4<-data%>%filter(target_name!="France")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Target = "Drop France") 

fance_tar5<-data%>%filter(target_name!="France")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Target = "Drop France") 

fance_tar6<-data%>%filter(target_name!="France")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Target = "Drop France") 


#Japan
japan_tar1<-data%>%filter(target_name!="Japan")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Target = "Drop Japan") 

japan_tar2<-data%>%filter(target_name!="Japan")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Target = "Drop Japan") 

japan_tar3<-data%>%filter(target_name!="Japan")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Target = "Drop Japan") 

japan_tar4<-data%>%filter(target_name!="Japan")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Target = "Drop Japan") 

japan_tar5<-data%>%filter(target_name!="Japan")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Target = "Drop Japan") 

japan_tar6<-data%>%filter(target_name!="Japan")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Target = "Drop Japan") 


#Israel
israel_tar1<-data%>%filter(target_name!="Israel")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Target = "Drop Israel") 

israel_tar2<-data%>%filter(target_name!="Israel")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Target = "Drop Israel") 

israel_tar3<-data%>%filter(target_name!="Israel")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Target = "Drop Israel") 

israel_tar4<-data%>%filter(target_name!="Israel")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Target = "Drop Israel") 

israel_tar5<-data%>%filter(target_name!="Israel")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Target = "Drop Israel") 

israel_tar6<-data%>%filter(target_name!="Israel")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Target = "Drop Israel") 

#China
china_tar1<-data%>%filter(target_name!="China")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Target = "Drop China") 

china_tar2<-data%>%filter(target_name!="China")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Target = "Drop China") 

china_tar3<-data%>%filter(target_name!="China")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Target = "Drop China") 

china_tar4<-data%>%filter(target_name!="China")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Target = "Drop China") 

china_tar5<-data%>%filter(target_name!="China")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Target = "Drop China") 

china_tar6<-data%>%filter(target_name!="China")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Target = "Drop China") 

#Iran
iran_tar1<-data%>%filter(target_name!="Iran")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Target = "Drop Iran") 

iran_tar2<-data%>%filter(target_name!="Iran")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Target = "Drop Iran") 

iran_tar3<-data%>%filter(target_name!="Iran")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Target = "Drop Iran") 

iran_tar4<-data%>%filter(target_name!="Iran")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Target = "Drop Iran") 

iran_tar5<-data%>%filter(target_name!="Iran")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Target = "Drop Iran") 

iran_tar6<-data%>%filter(target_name!="Iran")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Target = "Drop Iran") 

#North Korea
NK_tar1<-data%>%filter(target_name!="North Korea")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Target = "Drop North Korea") 

NK_tar2<-data%>%filter(target_name!="North Korea")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Target = "Drop North Korea") 

NK_tar3<-data%>%filter(target_name!="North Korea")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Target = "Drop North Korea") 

NK_tar4<-data%>%filter(target_name!="North Korea")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Target = "Drop North Korea") 

NK_tar5<-data%>%filter(target_name!="North Korea")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Target = "Drop North Korea") 

NK_tar6<-data%>%filter(target_name!="North Korea")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Target = "Drop North Korea") 

#Russia
russia_tar1<-data%>%filter(target_name!="Russia")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Target = "Drop Russia") 

russia_tar2<-data%>%filter(target_name!="Russia")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Target = "Drop Russia") 

russia_tar3<-data%>%filter(target_name!="Russia")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Target = "Drop Russia") 

russia_tar4<-data%>%filter(target_name!="Russia")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Target = "Drop Russia") 

russia_tar5<-data%>%filter(target_name!="Russia")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Target = "Drop Russia") 

russia_tar6<-data%>%filter(target_name!="Russia")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Target = "Drop Russia") 

#Venezuela
venezuela_tar1<-data%>%filter(target_name!="Venezuela")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Target = "Drop Venezuela") 

venezuela_tar2<-data%>%filter(target_name!="Venezuela")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Target = "Drop Venezuela") 

venezuela_tar3<-data%>%filter(target_name!="Venezuela")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Target = "Drop Venezuela") 

venezuela_tar4<-data%>%filter(target_name!="Venezuela")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Target = "Drop Venezuela") 

venezuela_tar5<-data%>%filter(target_name!="Venezuela")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Target = "Drop Venezuela") 

venezuela_tar6<-data%>%filter(target_name!="Venezuela")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Target = "Drop Venezuela") 

order_x <- c("Drop Australia", "Drop Canada", "Drop France", "Drop Japan", 
             "Drop Israel", "Drop China", "Drop Iran", 
             "Drop North Korea", "Drop Russia", "Drop Venezuela") 


#plot
dropcountry_target <- rbind(australia_tar1, australia_tar2,australia_tar3,
                          australia_tar4, australia_tar5,australia_tar6,
                          canada_tar1, canada_tar2, canada_tar3,
                          canada_tar4, canada_tar5, canada_tar6,
                          fance_tar1, fance_tar2, fance_tar3,
                          fance_tar4, fance_tar5, fance_tar6,
                          japan_tar1, japan_tar2, japan_tar3,
                          japan_tar4, japan_tar5, japan_tar6,
                          israel_tar1, israel_tar2, israel_tar3,
                          israel_tar4, israel_tar5, israel_tar6,
                          china_tar1, china_tar2, china_tar3,
                          china_tar4, china_tar5, china_tar6,
                          iran_tar1, iran_tar2, iran_tar3,
                          iran_tar4, iran_tar5, iran_tar6,
                          NK_tar1, NK_tar2, NK_tar3, 
                          NK_tar4, NK_tar5, NK_tar6, 
                          russia_tar1, russia_tar2, russia_tar3,
                          russia_tar4, russia_tar5, russia_tar6,
                          venezuela_tar1, venezuela_tar2, venezuela_tar3,
                          venezuela_tar4, venezuela_tar5, venezuela_tar6) %>% 
  filter(.,
         term == "shaming")

dropcountry_target$Outcome <- factor(dropcountry_target$Outcome, levels = c("Target Favorability", 
                                                                        "Military Cooperation w. Target", 
                                                                        "Economic Cooperation w. Target",
                                                                        "Shamer Favorability",
                                                                        "Military Cooperation w. Shamer", 
                                                                        "Economic Cooperation w. Shamer"))

# Generate coefficient plot 
ggplot(dropcountry_target, aes(x = Target, y = estimate, color = Outcome, shape = Outcome)) +
  geom_hline(yintercept = 0, color = "gray50", linetype = 2, size = 0.2) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.4)) +
  scale_colour_grey()+
  #coord_flip()+
  labs(x = "Model by dropped target state",
       y = "Effect Size of Shaming") +
  scale_x_discrete(limits= order_x) +
  theme(text = element_text(size = 20, family = "Times"),
        legend.key=element_blank(),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(angle = 45, hjust = 1, size = 15),
        plot.caption = element_text(size = 20, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

ggsave("Figures/dropcountry_target.pdf", width = 15, height = 5)


#Plot main effect of shaming by *shamer* country

#Australia
australia_sham1<-data%>%filter(shamer_name!="Australia")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Shamer = "Drop Australia") 

australia_sham2<-data%>%filter(shamer_name!="Australia")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Shamer = "Drop Australia") 

australia_sham3<-data%>%filter(shamer_name!="Australia")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Shamer = "Drop Australia") 

australia_sham4<-data%>%filter(shamer_name!="Australia")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Shamer = "Drop Australia") 

australia_sham5<-data%>%filter(shamer_name!="Australia")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Shamer = "Drop Australia") 

australia_sham6<-data%>%filter(shamer_name!="Australia")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Shamer = "Drop Australia") 

#Canada
canada_sham1<-data%>%filter(shamer_name!="Canada")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Shamer = "Drop Canada") 

canada_sham2<-data%>%filter(shamer_name!="Canada")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Shamer = "Drop Canada") 

canada_sham3<-data%>%filter(shamer_name!="Canada")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Shamer = "Drop Canada") 


canada_sham4<-data%>%filter(shamer_name!="Canada")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Shamer = "Drop Canada") 

canada_sham5<-data%>%filter(shamer_name!="Canada")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Shamer = "Drop Canada") 

canada_sham6<-data%>%filter(shamer_name!="Canada")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Shamer = "Drop Canada") 


#France
fance_sham1<-data%>%filter(shamer_name!="France")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Shamer = "Drop France") 

fance_sham2<-data%>%filter(shamer_name!="France")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Shamer = "Drop France") 

fance_sham3<-data%>%filter(shamer_name!="France")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Shamer = "Drop France") 

fance_sham4<-data%>%filter(shamer_name!="France")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Shamer = "Drop France") 

fance_sham5<-data%>%filter(shamer_name!="France")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Shamer = "Drop France") 

fance_sham6<-data%>%filter(shamer_name!="France")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Shamer = "Drop France") 


#Japan
japan_sham1<-data%>%filter(shamer_name!="Japan")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Shamer = "Drop Japan") 

japan_sham2<-data%>%filter(shamer_name!="Japan")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Shamer = "Drop Japan") 

japan_sham3<-data%>%filter(shamer_name!="Japan")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Shamer = "Drop Japan") 

japan_sham4<-data%>%filter(shamer_name!="Japan")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Shamer = "Drop Japan") 

japan_sham5<-data%>%filter(shamer_name!="Japan")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Shamer = "Drop Japan") 

japan_sham6<-data%>%filter(shamer_name!="Japan")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Shamer = "Drop Japan") 


#Israel
israel_sham1<-data%>%filter(shamer_name!="Israel")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Shamer = "Drop Israel") 

israel_sham2<-data%>%filter(shamer_name!="Israel")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Shamer = "Drop Israel") 

israel_sham3<-data%>%filter(shamer_name!="Israel")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Shamer = "Drop Israel") 

israel_sham4<-data%>%filter(shamer_name!="Israel")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Shamer = "Drop Israel") 

israel_sham5<-data%>%filter(shamer_name!="Israel")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Shamer = "Drop Israel") 

israel_sham6<-data%>%filter(shamer_name!="Israel")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Shamer = "Drop Israel") 

#China
china_sham1<-data%>%filter(shamer_name!="China")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Shamer = "Drop China") 

china_sham2<-data%>%filter(shamer_name!="China")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Shamer = "Drop China") 

china_sham3<-data%>%filter(shamer_name!="China")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Shamer = "Drop China") 

china_sham4<-data%>%filter(shamer_name!="China")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Shamer = "Drop China") 

china_sham5<-data%>%filter(shamer_name!="China")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Shamer = "Drop China") 

china_sham6<-data%>%filter(shamer_name!="China")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Shamer = "Drop China") 

#Iran
iran_sham1<-data%>%filter(shamer_name!="Iran")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Shamer = "Drop Iran") 

iran_sham2<-data%>%filter(shamer_name!="Iran")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Shamer = "Drop Iran") 

iran_sham3<-data%>%filter(shamer_name!="Iran")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Shamer = "Drop Iran") 

iran_sham4<-data%>%filter(shamer_name!="Iran")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Shamer = "Drop Iran") 

iran_sham5<-data%>%filter(shamer_name!="Iran")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Shamer = "Drop Iran") 

iran_sham6<-data%>%filter(shamer_name!="Iran")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Shamer = "Drop Iran") 

#North Korea
NK_sham1<-data%>%filter(shamer_name!="North Korea")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Shamer = "Drop North Korea") 

NK_sham2<-data%>%filter(shamer_name!="North Korea")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Shamer = "Drop North Korea") 

NK_sham3<-data%>%filter(shamer_name!="North Korea")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Shamer = "Drop North Korea") 

NK_sham4<-data%>%filter(shamer_name!="North Korea")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Shamer = "Drop North Korea") 

NK_sham5<-data%>%filter(shamer_name!="North Korea")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Shamer = "Drop North Korea") 

NK_sham6<-data%>%filter(shamer_name!="North Korea")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Shamer = "Drop North Korea") 

#Russia
russia_sham1<-data%>%filter(shamer_name!="Russia")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Shamer = "Drop Russia") 

russia_sham2<-data%>%filter(shamer_name!="Russia")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Shamer = "Drop Russia") 

russia_sham3<-data%>%filter(shamer_name!="Russia")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Shamer = "Drop Russia") 

russia_sham4<-data%>%filter(shamer_name!="Russia")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Shamer = "Drop Russia") 

russia_sham5<-data%>%filter(shamer_name!="Russia")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Shamer = "Drop Russia") 

russia_sham6<-data%>%filter(shamer_name!="Russia")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Shamer = "Drop Russia") 

#Venezuela
venezuela_sham1<-data%>%filter(shamer_name!="Venezuela")%>%
  lm_robust(target_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Target Favorability",
         Shamer = "Drop Venezuela") 

venezuela_sham2<-data%>%filter(shamer_name!="Venezuela")%>%
  lm_robust(target_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Target",
         Shamer = "Drop Venezuela") 

venezuela_sham3<-data%>%filter(shamer_name!="Venezuela")%>%
  lm_robust(target_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Target",
         Shamer = "Drop Venezuela") 

venezuela_sham4<-data%>%filter(shamer_name!="Venezuela")%>%
  lm_robust(shamer_favorability~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Shamer Favorability",
         Shamer = "Drop Venezuela") 

venezuela_sham5<-data%>%filter(shamer_name!="Venezuela")%>%
  lm_robust(shamer_military~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Military Cooperation w. Shamer",
         Shamer = "Drop Venezuela") 

venezuela_sham6<-data%>%filter(shamer_name!="Venezuela")%>%
  lm_robust(shamer_economic~shaming, data = .)%>% tidy %>% 
  mutate(.,
         Outcome = "Economic Cooperation w. Shamer",
         Shamer = "Drop Venezuela") 

order_x <- c("Drop Australia", "Drop Canada", "Drop France", "Drop Japan", 
             "Drop Israel", "Drop China", "Drop Iran", 
             "Drop North Korea", "Drop Russia", "Drop Venezuela") 


#plot
dropcountry_shamer <- rbind(australia_sham1, australia_sham2,australia_sham3,
                          australia_sham4, australia_sham5,australia_sham6,
                          canada_sham1, canada_sham2, canada_sham3,
                          canada_sham4, canada_sham5, canada_sham6,
                          fance_sham1, fance_sham2, fance_sham3,
                          fance_sham4, fance_sham5, fance_sham6,
                          japan_sham1, japan_sham2, japan_sham3,
                          japan_sham4, japan_sham5, japan_sham6,
                          israel_sham1, israel_sham2, israel_sham3,
                          israel_sham4, israel_sham5, israel_sham6,
                          china_sham1, china_sham2, china_sham3,
                          china_sham4, china_sham5, china_sham6,
                          iran_sham1, iran_sham2, iran_sham3,
                          iran_sham4, iran_sham5, iran_sham6,
                          NK_sham1, NK_sham2, NK_sham3, 
                          NK_sham4, NK_sham5, NK_sham6, 
                          russia_sham1, russia_sham2, russia_sham3,
                          russia_sham4, russia_sham5, russia_sham6,
                          venezuela_sham1, venezuela_sham2, venezuela_sham3,
                          venezuela_sham4, venezuela_sham5, venezuela_sham6) %>% 
  filter(.,
         term == "shaming")

dropcountry_shamer$Outcome <- factor(dropcountry_shamer$Outcome, levels = c("Target Favorability", 
                                                                        "Military Cooperation w. Target", 
                                                                        "Economic Cooperation w. Target",
                                                                        "Shamer Favorability",
                                                                        "Military Cooperation w. Shamer", 
                                                                        "Economic Cooperation w. Shamer"))

# Generate coefficient plot 
ggplot(dropcountry_shamer, aes(x = Shamer, y = estimate, color = Outcome, shape = Outcome)) +
  geom_hline(yintercept = 0, color = "gray50", linetype = 2, size = 0.2) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.4)) +
  scale_colour_grey()+
  #coord_flip()+
  labs(x = "Model by dropped shamer state",
       y = "Effect Size of Shaming") +
  scale_x_discrete(limits= order_x) +
  theme(text = element_text(size = 20, family = "Times"),
        legend.key=element_blank(),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(angle = 45, hjust = 1, size = 15),
        plot.caption = element_text(size = 20, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

ggsave("Figures/dropcountry_shamer.pdf", width = 15, height = 5)



#probing external validity bias using Devaux and Egami
#report for main significant findings
covariates <- c("Democrat","Republican","Independent", "Educated", 
                "Black","White", "Asian", "Native", "Other", "Religious", "Female", 
                "Age", "Political_interest", "Conservative", "hawk1", "hawk2", "hawk3", 
                "hawk4","state") 

set.seed(100)
#effect of shaming on target outcomes
#targer_favorability
for_sate1 <- as.formula(paste0("target_favorability ~ shaming + ", paste(covariates, collapse = "+")))
lm_sate1 <- lm(for_sate1, data = data)
sate_est1 <- summary(lm_sate1)$coef["shaming", c(1,2)]
sate_est1

exr_out1 <- exr(outcome = "target_favorability", 
                treatment = "shaming", 
                covariates = covariates, 
                data = data,
                boot = TRUE,
                sate_estimate = sate_est1)


#target_military
for_sate1a <- as.formula(paste0("target_military ~ shaming + ", paste(covariates, collapse = "+")))
lm_sate1a <- lm(for_sate1a, data = data)
sate_est1a <- summary(lm_sate1a)$coef["shaming", c(1,2)]
sate_est1a

exr_out1a <- exr(outcome = "target_military", 
                 treatment = "shaming", 
                 covariates = covariates, 
                 data = data,
                 boot = TRUE,
                 sate_estimate = sate_est1a) 

#target_economy
for_sate1b <- as.formula(paste0("target_economic ~ shaming + ", paste(covariates, collapse = "+")))
lm_sate1b <- lm(for_sate1b, data = data)
sate_est1b <- summary(lm_sate1b)$coef["shaming", c(1,2)]
sate_est1b

exr_out1b <- exr(outcome = "target_economic", 
                 treatment = "shaming", 
                 covariates = covariates, 
                 data = data,
                 boot = TRUE,
                 sate_estimate = sate_est1b) 


#effect of shaming on shamer outcomes
#shamer_fav
for_sate2 <- as.formula(paste0("shamer_favorability ~ shaming + ", paste(covariates, collapse = "+")))
lm_sate2 <- lm(for_sate2, data = data)
sate_est2 <- summary(lm_sate2)$coef["shaming", c(1,2)]
sate_est2

exr_out2 <- exr(outcome = "shamer_favorability", 
                treatment = "shaming", 
                covariates = covariates, 
                data = data,
                boot = TRUE,
                sate_estimate = sate_est2)
#military coop
for_sate2a <- as.formula(paste0("shamer_military ~ shaming + ", paste(covariates, collapse = "+")))
lm_sate2a <- lm(for_sate2a, data = data)
sate_est2a <- summary(lm_sate2a)$coef["shaming", c(1,2)]
sate_est2a

exr_out2a <- exr(outcome = "shamer_military", 
                 treatment = "shaming", 
                 covariates = covariates, 
                 data = data,
                 boot = TRUE,
                 sate_estimate = sate_est2a)

#economic coop
for_sate2b <- as.formula(paste0("shamer_economic ~ shaming + ", paste(covariates, collapse = "+")))
lm_sate2b <- lm(for_sate2b, data = data)
sate_est2b <- summary(lm_sate2b)$coef["shaming", c(1,2)]
sate_est2b

exr_out2b <- exr(outcome = "shamer_economic", 
                 treatment = "shaming", 
                 covariates = covariates, 
                 data = data,
                 boot = TRUE,
                 sate_estimate = sate_est2b)


#arbitrary detention
for_sate3a <- as.formula(paste0("detention_new ~ shaming + ", paste(covariates, collapse = "+")))
lm_sate3a <- lm(for_sate3a, data = data)
sate_est3a <- summary(lm_sate3a)$coef["shaming", c(1,2)]
sate_est3a

exr_out3a <- exr(outcome = "detention_new", 
                 treatment = "shaming", 
                 covariates = covariates, 
                 data = data,
                 boot = TRUE,
                 sate_estimate = sate_est3a)


#support HRs
for_sate3b <- as.formula(paste0("HRs_new ~ shaming + ", paste(covariates, collapse = "+")))
lm_sate3b <- lm(for_sate3b, data = data)
sate_est3b <- summary(lm_sate3b)$coef["shaming", c(1,2)]
sate_est3b

exr_out3b <- exr(outcome = "HRs_new", 
                 treatment = "shaming", 
                 covariates = covariates, 
                 data = data,
                 boot = TRUE,
                 sate_estimate = sate_est3b)

#index
for_sate3c <- as.formula(paste0("H3_index ~ shaming + ", paste(covariates, collapse = "+")))
lm_sate3c <- lm(for_sate3c, data = data)
sate_est3c <- summary(lm_sate3c)$coef["shaming", c(1,2)]
sate_est3c

exr_out3c <- exr(outcome = "H3_index", 
                 treatment = "shaming", 
                 covariates = covariates, 
                 data = data,
                 boot = TRUE,
                 sate_estimate = sate_est3c)



pdf(file="Figures/extr_general.pdf", width = 15, height = 10)
par(mfrow=c(3,3))

plot(exr_out1)
mtext("Target favoraility",side=3)
plot(exr_out1a)
mtext("Target military coop.",side=3)
plot(exr_out1b)
mtext("Target economic coop.",side=3)

plot(exr_out2)
mtext("Shamer favorability",side=3)
plot(exr_out2a)
mtext("Shamer military coop.",side=3)
plot(exr_out2b)
mtext("Shamer economic coop.",side=3)

plot(exr_out3a)
mtext("Oppose arbitrary detention",side=3)
plot(exr_out3b)
mtext("Support HRs",side=3)
plot(exr_out3c)
mtext("H3 Index",side=3)
dev.off()



#Plot outcome by condition


plot_outcome <- function(data, outcome_var, y_label) {
  
  # Prepare data
  data_new <- data %>%
    mutate(
      !!outcome_var := as.numeric(.data[[outcome_var]]),
      condition_2 = case_when(
        shaming == 1 ~ "Shaming",
        shaming == 0 ~ "Control"
      )
    ) %>%
    filter(
      !is.na(condition_2),
      !is.na(.data[[outcome_var]]),
      .data[[outcome_var]] >= 1,
      .data[[outcome_var]] <= 5
    )
  
  # Group summaries
  summary_df <- data_new %>%
    group_by(condition_2) %>%
    do(tidy(lm_robust(reformulate("1", outcome_var), data = .))) %>%
    mutate(value = estimate) %>%
    drop_na()
  
  # p-value
  model <- lm_robust(reformulate("condition_2", outcome_var), data = data_new)
  p_val <- model$p.value[2]
  p_label <- if (p_val < 0.01) {
    "italic(p) < 0.01"
  } else {
    paste0("italic(p) == ", format(round(p_val, 3), nsmall = 3))
  }
  
  # Plot
  ggplot(summary_df, aes(condition_2, value)) +
    # Smaller, more dispersed jittered dots
    geom_point(
      data = data_new,
      aes(x = condition_2, y = .data[[outcome_var]]),
      position = position_jitter(width = 0.2, height = 0.3),
      alpha = 0.2,
      size = 0.7
    ) +
    # Means and error bars, nudged right
    geom_point(size = 3, color = "firebrick3", position = position_nudge(x = 0.1)) +
    geom_line(color = "firebrick3", linetype = "dashed", group = 1, position = position_nudge(x = 0.1)) +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, color = "firebrick3", position = position_nudge(x = 0.1)) +
    geom_text(aes(label = round(estimate, 2)), vjust = 1.2, hjust = -1, size = 4, color = "firebrick3", position = position_nudge(x = 0.1)) +
    # P-value
    annotate("text", x = 2.4, y = 1.1, label = p_label, parse = TRUE, family = "Times", size = 4, color = "darkblue") +
    coord_cartesian(ylim = c(1, 5)) +
    theme_bw() +
    ylab(y_label) +
    theme(axis.title.x = element_blank())
}


#target outcomes by condition
p1 <- plot_outcome(data, "target_favorability", "Target Favorability")
p2 <- plot_outcome(data, "target_military", "Support for Military Cooperation with Target")
p3 <- plot_outcome(data, "target_economic", "Support for Economic Cooperation with Target")

pdf("Figures/target_bycond.pdf", width = 12, height = 4)
(p1 | p2 | p3) + plot_layout(guides = "collect")
dev.off()


#shamer outcomes by condition
p4 <- plot_outcome(data, "shamer_favorability", "Shamer Favorability")
p5 <- plot_outcome(data, "shamer_military", "Support for Military Cooperation with Shamer")
p6 <- plot_outcome(data, "shamer_economic", "Support for Economic Cooperation with Shamer")

pdf("Figures/shamer_bycond.pdf", width = 12, height = 4)
(p4 | p5 | p6) + plot_layout(guides = "collect")
dev.off()


#HRs outcomes by condition
p7 <- plot_outcome(data, "detention_new", "Oppose Arbitrary Detention")
p8 <- plot_outcome(data, "HRs_new", "Support for Human Rights")
p9 <- plot_outcome(data, "law_new", "Support International Law")
p10 <- plot_outcome(data, "H3_index", "Index")


pdf("Figures/HRs_bycond.pdf", width = 16, height = 4)
(p7 | p8 | p9| p10) + plot_layout(guides = "collect")
dev.off()



